# ===============================================
# Paris ratchet paper
# Code to run extensions of main models
# ===============================================

## ## Contents

#' 1. Load data for models from "00_make_analysis_data.R"
#' 2. Regression on timing of submissions
#' 3. Measurement and accounting issues
#' 4. Climate laws

#### Initialize #### 

library(tidyverse)
library(knitr)
library(DescTools) # Winsorize()
library(broom)
library(modelsummary)
library(nnet) # for multinomial logit
library(magrittr)
library(lubridate)


## Set working directory


# Define core ISO3C codes to use in analysis
core_iso3c <- c("AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEL", "BEN", "BFA", "BGD", "BGR", "BHR", "BHS", "BIH", "BLR", "BLZ", "BOL", "BRA", "BRB", "BRN", "BTN", "BWA", "CAF", "CAN", "CHE", "CHL", "CHN", "CIV", "CMR", "COD", "COG", "COK", "COL", "COM", "CPV", "CRI", "CUB", "CYP", "CZE", "DEU", "DJI", "DMA", "DNK", "DOM", "DZA", "ECU", "EGY", "ERI", "ESP", "EST", "ETH", "FIN", "FJI", "FRA", "FSM", "GAB", "GBR", "GEO", "GHA", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GTM", "GUY", "HND", "HRV", "HTI", "HUN", "IDN", "IND", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAO", "LBN", "LBR", "LBY", "LCA", "LIE", "LKA", "LSO", "LTU", "LUX", "LVA", "MAR", "MDA", "MDG", "MDV", "MEX", "MHL", "MKD", "MLI", "MLT", "MMR", "MNE", "MNG", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAM", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRU", "NZL", "OMN", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRK", "PRT", "PRY", "QAT", "ROU", "RUS", "RWA", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SOM", "SRB", "SSD", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SYC", "SYR", "TCD", "TGO", "THA", "TJK", "TKM", "TLS", "TON", "TTO", "TUN", "TUR", "TUV", "TZA", "UGA", "UKR", "URY", "USA", "UZB", "VCT", "VEN", "VNM", "VUT", "WSM", "YEM", "ZAF", "ZMB", "ZWE")

## Load data
#' code to make these objects in "00_make_analysis_data.R"

country_level_average <- readRDS("output/country_level_average_20230920.rds")

data_for_models <- readRDS("output/data_for_models20230920.Rds")




#### Regressions on timing of submission ####

eu_iso <- c("AUT", "BEL", "BGR", "HRV", "CYP", 
            "CZE", "DNK", "EST", "FIN", "FRA", 
            "DEU", "GRC", "HUN", "IRL", "ITA", 
            "LVA", "LTU", "LUX", "MLT", "NLD", 
            "POL", "PRT", "ROU", "SVK", "SVN", 
            "ESP", "SWE")

data_for_date <- country_level_average %>% 
  select(iso3c, date = ndc2_date) %>% 
  distinct() %>% 
  mutate(date = replace(date, 
                        iso3c %in% eu_iso, 
                        as_date("2020-12-18"))) %>% 
  mutate(date = replace(date, 
                        is.na(date), 
                        as_date("2022-09-01"))) %>% 
  mutate(update_ndc_buckets = case_when(date == as_date("2022-09-01") ~ 0, # Countries that never submit
                                        date < as_date("2019-01-01") ~ 0, # Countries that submit way too early
                                        date > as_date("2019-01-01") & date <= as_date("2021-01-01") ~ 1, # Countries that submit pre-EU
                                        date > as_date("2021-01-01") & date <= as_date("2021-05-01") ~ 2, # Countries that submit post-EU, pre-US
                                        date > as_date("2021-05-01") ~ 3)) %>% # Countries that submit post-US
  mutate(update_ndc_labels = as.factor(case_when(update_ndc_buckets == 0 ~ "Never",
                                                 update_ndc_buckets == 1 ~ "Early",
                                                 update_ndc_buckets == 2 ~ "Post-EU",
                                                 update_ndc_buckets == 3 ~ "Post-US"))) %>% 
  # mutate(update_ndc_labels = relevel(update_ndc_labels, ref = "Never")) %$% table(update_ndc_labels)
  mutate(update_ndc_labels = relevel(update_ndc_labels, ref = "Never")) %>%
  full_join(., data_for_models, by = "iso3c") %>% 
  filter(!(iso3c %in% eu_iso), 
         !(iso3c %in% c("USA", "NOR", "ISL", "CHE"))) 

multinomial_logit_date <- data_for_date %>% 
  multinom(update_ndc_labels ~ 
             trade_exposure_USA + 
             trade_exposure_EU, 
           data = .) 


table_mlogit_timing <- multinomial_logit_date %>% 
  tidy() %>% 
  mutate(order = case_when(y.level == "Early" ~ 1,
                           y.level == "Post-EU" ~ 2,
                           y.level == "Post-US" ~ 3)) %>% 
  arrange(order, term) %>% 
  mutate_if(is.numeric, ~round(., 3)) %>% 
  mutate(estimate = ifelse(p.value < 0.05, 
                           paste(estimate, "*", sep = ""), estimate)) %>% 
  select(-order, -statistic, -p.value) %>% 
  as.data.frame() %>% 
  mutate(term = case_when(term == "trade_exposure_EU" ~ "EU trade exposure",
                          term == "trade_exposure_CHN" ~ "China trade exposure",
                          term == "trade_exposure_USA" ~ "US trade exposure",
                          term == "(Intercept)" ~ "(Intercept)")) %>% 
  pivot_wider(names_from = "y.level", values_from = "estimate":"std.error") %>% 
  set_colnames(c("term", "early_beta", "post_eu_beta", "post_us_beta", 
                 "early_stderror", "post_eu_stderror", "post_us_stderror")) %>% 
  select(1, 2, 5, 3, 6, 4, 7) %>% 
  knitr::kable(., format = "latex", booktabs = T)
table_mlogit_timing
#' -> paste into manuscript
#' Check caption


#### Measurement and accounting issues ####

## Countries with missing targets

m_imputed1 <- data_for_models %>% 
  lm(glasgow_target_impute ~ paris_target_impute + 
       io_percentage + tradeflow_percentage, 
     data = .)


## Models without imputed targets in building spatial weights
m_tradeflow_safe <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe + 
       io_safe_percentage + tradeflow_safe_percentage, 
     data = .)

## Countries with outlying target 
m_outliers1 <- data_for_models %>% 
  filter(!(iso3c %in% c("GMB", "COD", "PRY"))) %>%
  lm(glasgow_target_safe ~ paris_target_safe + 
       io_percentage + tradeflow_percentage, 
     data = .)

## Outliers in main model based on Cook's distance
m1 <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage, 
     data = .)
# cooks.distance(m1)
to_drop_cooks <- cbind.data.frame(
  leverage = 
    cooks.distance(m1) > (4/length(cooks.distance(m1)))) %>% 
  filter(leverage == T) %>% 
  rownames_to_column() %>% 
  pull(1)

data_for_models %>% 
  slice(as.integer(to_drop_cooks)) %>% 
  select(iso3c)

m_outliers2 <- data_for_models %>%
  filter(!row_number() %in% to_drop_cooks) %>% 
  lm(glasgow_target_safe ~ paris_target_safe + 
       io_percentage + tradeflow_percentage, 
     data = .)

## Models with PRIMAP GHG series

# Baseline model, to check for outliers
m_primap_full <- data_for_models %>% 
  lm(glasgow_target_primap ~ paris_target_primap +
       io_percentage + tradeflow_percentage,
     data = .)

# Identify influential outliers
to_drop_cooks_primap <- cbind.data.frame(leverage = cooks.distance(m_primap_full) > 
                                           (4/length(cooks.distance(m_primap_full)))) %>% 
  filter(leverage == T) %>% 
  rownames_to_column() %>% 
  pull(1)

m_primap_cooks <- data_for_models %>% 
  filter(!row_number() %in% to_drop_cooks_primap) %>% 
  # filter(!iso3c %in% c("ARM", "BFA", "BTN", "COD", "GEO", "GMB", "KHM", "KNA", "MLI", "VNM")) %>%
  lm(glasgow_target_primap ~ paris_target_primap +
       io_percentage + tradeflow_percentage,
     data = .)

gm <- tibble::tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Observations", 0,
  "r.squared", "R2", 2)
coefficients_renamed = c(
  "tradeflow_percentage" = "Trade Paris",
  "eite_percentage" = "Trade (EITE) Paris",
  "green_percentage" = "Trade (Green) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "io_percentage" = "IO Paris",
  "io_safe_percentage" = "IO Paris (don't impute partner's target)",
  "tradeflow_safe_percentage" = "Trade Paris (don't impute partner's target)",
  "paris_target_primap" = "Paris target (Excl. LULUCF)",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "paris_target_safe" = "Paris target",
  "paris_target_impute" = "Paris target (impute)",
  "(Intercept)" = "(Intercept)")

table_targets_measurement <- modelsummary(
  models = list(m_imputed1, m_tradeflow_safe, 
                m_outliers1, m_outliers2, m_primap_cooks),
             output = "latex",
             fmt = 2,
             coef_rename = coefficients_renamed,
             gof_map = gm, 
             notes = list('Outcome is Glasgow mitigation target'),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table_targets_measurement
#' -> paste into manuscript
#' Add "\midrule" before controls 
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption



# #### Climate laws ####

m_laws1 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       tradeflow_laws,
     data = .)

m_laws2 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       tradeflow_laws_new,
     data = .)

m_laws3 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       io_laws,
     data = .)

m_laws4 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       io_laws_new,
     data = .)

m_laws5 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       tradeflow_laws + io_laws,
     data = .)

m_laws6 <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe +
       # trade_openness + industry + renewables_share + fossil_rents +
       tradeflow_laws_new + io_laws_new,
     data = .)

gm <- tibble::tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Observations", 0,
  "r.squared", "R2", 2)

coefficients_renamed = c(
  "tradeflow_percentage" = "Trade Paris",
  "eite_percentage" = "Trade (EITE) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "io_percentage" = "IO Paris",
  "io_safe_percentage" = "IO Paris (don't impute partner's target)",
  "tradeflow_safe_percentage" = "Trade Paris (don't impute partner's target)",
  "paris_target_primap" = "Paris target (Excl. LULUCF)",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "paris_target_safe" = "Paris target",
  "paris_target_impute" = "Paris target (impute)",
  #
  "tradeflow_laws" = "Trade Laws (count)",
  "tradeflow_laws_new" = "Trade Laws (post-Paris)",
  "io_laws" = "IO Laws (count)",
  "io_laws_new" = "IO Laws (post-Paris)",
  "(Intercept)" = "(Intercept)")


table_targets_laws <- modelsummary(
  models = list(m_laws1, m_laws2, m_laws3, m_laws4,
                m_laws5, m_laws6),
             output = "latex",
             coef_rename = coefficients_renamed,
             gof_map = gm,
             fmt = 2,
             notes = list('Outcome is Glasgow mitigation target'),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))

table_targets_laws
#' -> paste into manuscript
#' Add "\midrule" before controls 
#' Bold the spatial matrix terms
#' Math the R2 term
#' Check caption




